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SECTION  I 
INTRODUCTION 


The  problem  of  inferring  a  mechanical  relaxation  spectrum 
in  linear  viscoelastic  theory  has  been  the  subject  of  many  papers 
(References  1-8) .  The  experimentally  obtained  modulus  functions 
are  related  to  the  relaxation  spectriun  through  the  equations 
(Reference  9) 


II 

1  H(t) 

exp(-t/T)d  log  t. 
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1) 

E'(W)  =  J 

^  00 

f  H(t) 
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3) 

where  E^(t)  is  the  stress  relaxation  modulus,  E'(w)  and  E”(w) 
are  the  in-phase  and  out-of-phase  components  of  the  complex 
dynamic  modulus,  respectively,  the  H(t)  is  the  sought-after 
relaxation  spectrum. 


Any  one  of  the  above  equations  can  be  expressed  in  the 
general  form  of 


'/^(x)  = 


■f 


K(x,y)  f  (y)dy 


(  4) 


or,  in  matrix  notation 

^  =  g  f  (  5) 

where  ^  is  the  measured  experimental  data  of  the  problem 
resulting  in  a  solution  f  and  K  is  the  operator  (the  kernel 
in  this  case)  that  relates  the  two  vectors.  If  f®  denotes 
the  real  solution,  there  should  exist  an  "exact"  set  of 
{  t|j°(x)}  such  that 

i"  =  g  (6) 

Experimentally,  only  an  approximation  ^  to  can  be  obtained. 
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The  task  at  hand  is  then  to  find  a  procedure  to  infer  a 
solution  f  from  such  that  for  a  given  f  is  the  best 
approximation  of  f®.  This  means  we  must  find  an  appropriate 
operator  (or  procedure)  B  (not  necessary  ,  such  that 

f®  =  B  i  .  (7) 

B 

One  approach  is  to  find  a  solution  f  which  satisfies  the 
condition  that  the  function 

N  =  (J  -  K  f®)2  (8) 

to  be  a  minimum. 

If  the  problem  is  an  ill-posed  problem,  the  solution  thus 
found  may  not  be  a  good  solution.  Tikhonov  (Reference  10)  has 
found  that  in  such  instances,  by  adding  a  regularization  term 
to  the  solution  mapping  procedure,  B,  a  different  f*  (the 
asterisk  denotes  the  solution  is  obtained  with  regularization) 
can  be  obtained  which  is  a  better  approximation  to  the  real 
solution  f ®  than  f  .  Using  a  scaling  factor  a ,  to  control  the 
magnitude  of  the  regularization  term  contribution  in  the  mapping 

procedure,  B,  it  was  shown  that  there  exists  an  optimum  a  which 

“  ~  2 
gives  a  minimiim  value  to  the  square  difference  of  (f*  -  f®)  . 

Wiff  and  Gehatia  (Reference  11)  have  combined  the  regulari¬ 
zation  method  and  a  quadratic  programming  algorithm  to  treat 
the  ill-posed  problem  of  determining  a  molecular  weight  distri¬ 
bution  of  polymers  from  sedimentation-diffusion  equilibrium  data 
obtained  via  ultracentrifuge  experiments.  Since  in  real  situations 
a  priori  knowledge  of  f®  is  not  available,  the  value  of  the  optimvim 
a  is  not  known.  However,  it  was  found  that  the  norm  | | (^  -  S  f*)  1 | 
will  go  through  a  minimum  with  increasing,  a.  This  solution 

f '  (the  prime  denoting  the  solution  that  gives  the  minimum 

'^B 

norm)  is  a  much  better  approximation  of  f®  than  f  .  Using  this 
approach  the  Regularization  Quadratic  Programming  (RQP)  method 
has  been  used  successfully  in  the  molecular  weight  distribution 
determination  (Reference  12)  . 
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Wiff  (Reference  1)  had  applied  the  RQP  method  to  the  problem 
of  inferring  relaxation  spectra  from  modulus  data  and  the  pre¬ 
liminary  results  were  promising.  This  work  is  a  more  in-depth 
study  of  the  problem.  It  is  found  that  regularization  will  give 
a  minimum  norm  of  | | ( J  -  K  i*)^| 1  only  if  different  integration 
routines  are  used  in  the  forward  mapping  operator  B  and  the  back 
calculation  operator  K.  And  because  of  the  nature  of  the  relax¬ 
ation  spectrum,  regularization  is  not  suitable  for  this  problem. 

A  modification  to  the  expression  to  be  minimized  (Equation  7) 
was  made.  By  adjusting  the  number  of  output  points,  a  reasonable 
approximation  of  a  known  relaxation  spectrum  has  been  obtained. 
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SECTION  II 

REGULARIZATION  QUADRATIC  PROGRAMMING  (RQP) 

Quadratic  programming  can  be  used  to  find  the  minimiam  of 
the  expression  N  (Equation  8) .  The  algorithm  (Reference  13) 
will  find  a  solution  x  that  minimizes  the  expression 

M  =  xgx-2Ex  (9) 

if  D  and  E  are  given.  Expanding  Equation  (9) ,  one  can  show 
(with  daggers  denoting  ad joints) 

N  =  f g"*"  g  f  -  2  J’*'  K  f  +  J  .  (10) 

Comparing  with  Equation  (9)  and  dropping  the  last  term  because 
it  is  a  constant,  it  is  obvious  that 


g  =  S  2  , 

^  ■'.i  ■'nj 


(Ha) 


E  =  i  K 

=  hi 

The  regularization  term  is 


fi(n) 


[f(y)l  = 


r 

/  = 

a  L 


d"f(v) 


(lib) 


By  minimizing  this  term  plus  the  functional  in  Equation  (10)  the 
solution  f*  is  smoothed  to  remove  noises  which  were  assiamed 
present  due  to  the  ill-posedness  of  the  problem.  Expressing  the 
derivative  in  terms  of  binomial  expansion: 


^n  \k/  ^  • 


^j-p+k 


where  h  is  the  const?nt  increment  in  y. 


0‘ 


is  the  binomial 


coefficient,  and  p=n  for  n=odd  and  p=n-l  for  n=even.  Substituting 
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Equation  (13)  in  the  matrix  form  of  Equation  (12) ,  the  regulari¬ 
zation  form  is  then 


[f] 


_1 _ 

h^-1 


j  j  n  n 

E  I  S  Z 

g=l  r=l  k=0  1=0 


(-1)^^^  f  f 
^  ’  g-p+k  r-p+1. 


(14) 


which  in  matrix  notation  will  be 

[f]  =  f  (15) 

The  RQP  method  then  uses  quadratic  programming  to  minimize 
the  function 

M  =  N  +  a^2  (16) 

where  N  and  are  defined  in  Equations  (10)  and  (15) ,  respect¬ 
ively,  and  a  is  a  scaling  factor  that  controls  the  relative  con¬ 
tributions  of  the  two  terms  to  the  total  function  M.  With  a  too 
small,  the  ill-posedness  of  M  was  postulated  to  result  in  too 
much  noise  in  the  solution  f*;  but  when  a  was  too  large,  the 
solution  became  over-smoothed  and  deviated  from  the  sought-after 
function  f ® . 

The  scheme  (I)  used  in  Reference  lU  varies  a  to  obtain  a 
set  of  f*,  then  the  norms 

I  1  (J  -  K  f*)^l  I  (17) 

were  calculated,  and  the  function  f '  that  gives  the  inferior  norm 
value  was  used  as  the  best  approximation  of  .  The  numerical 
integration  used  in  the  term  N  (Equation  16)  was  the  Trapezoidal 
Rule,  while  the  Simpson's  Rule  was  used  in  calculating  the  norm 
(Equation  17) . 

Modifications  to  this  program  have  been  made  during  this 
work.  The  new  scheme  (II)  has  two  major  differences  from  scheme 
(I) .  Firstly,  the  expression  N  is  changed  to  reflect  the  relative 
difference,  rather  than  the  linear  difference  of  the  data  points. 
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2 


The  new  expression  is 


i  i  s  i 


i  ^ 


s  i 

2  -  +  1 


(18) 


Comparing  with  Equation 


and 


(19a) 

(19b) 


Secondly,  the  integration  routine  in  the  calculation  of  norm  is 
changed  to  Trapezoidal  Rule,  to  be  consistent  with  the  routine 
used  in  Equation  (18) . 


In  both  schemes,  the  domain  of  the  solution  has  to  be 
specified.  The  range  of  t  for  the  relaxation  spectrum  should  be 
large  enough  to  cover  the  input  data  range.  If  t^^  and  t2  (u^ 
and  ^2  for  storage  and  loss  modulus  data)  are  the  lower  and  upper 
limit  of  the  input  data,  then  <  t^  and  ^2  ^  ^2 

and  >  1/^2^ •  This  is  because  the  kernel  for  each  value  of 
E(t)  covers  a  range  of  t.  For  E(tj^)  ,  contributions  from  H(t)  for 
T  <  t^  are  present.  By  setting  a  limit  of  the  solution  form 
is  forced  to  have  H(t)  =  0  for  t  <  if  no  H(t)  value  is 

allowed  for  ^1^  errors  will  be  forced  into  the  region  between 
and  to  account  for  the  relaxation  effect  in  E(tj^)  .  For  the 
same  reasoning,  the  H(t)  values  obtained  near  both  t  limits  should 
be  regarded  with  caution.  Unlike  other  experimental  data,  there 
is  no  guarantee  that  the  mechanical  modulus  data  is  complete,  and 
fully  account  for  the  changes  reflected  by  H(t) .  Often,  instru¬ 
ment  sensitivity  determines  the  range  of  the  input  data,  or  in  a 
master  curve  obtained  near  the  glass  transition  region,  the  limit 
in  the  short  time  region  is  where  the  experimenter  believes  that 
the  assinnptions  underlying  the  time-temperature  superposition 
principle  are  no  longer  valid. 
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SECTION  III 

RESULTS  AND  DISCUSSION 

The  published  experimental  data  of  Tobolsky  and  Catsiff 
[Reference  14]  on  poly isobutylene  were  used  to  check  the  capability 
of  scheme  (1)  to  determine  relaxation  spectrum  from  mechanical 
data.  Figure  l(a-c)  shows  the  H(t) ,  (f')  results  generated  from 
the  stress-relaxation,  storage  and  loss  moduli,  respectively,  and 
all  of  which  contain  unacceptable  levels  of  noise  in  the  solution 
forms . 

Figure  1 (d)  is  the  comparison  between  the  initial  stress- 

relaxation  data  and  the  back-calculated  modulus  values  using  the 

computer  generated  H(t) .  It  is  obvious  that  the  high  magnitude 

data  fit  better  than  those  of  lower  magnitude,  which  is  not 

surprising  because  the  initial  data  vary  in  magnitude  over  three 
10  7  2 

decades  (10  to  10  dynes/cm  ).  When  the  computer  algorithm 
minimizes  the  sum  in  Equation  (10),  the  terms  from  the  large 
magnitude  data  will  dominate  the  minimization  process.  By  changing 
the  sum  to  be  minimized  in  scheme  II  to  Equation  (18) ,  which  con¬ 
siders  the  relative  deviation  of  each  data  point  instead  of  the 
absolute  magnitude  of  deviations,  all  data  points  are  weighted 
equally  in  the  process. 

In  the  original  RQP  scheme,  without  regularization,  the 
quadratic  programming  algorithm  is  supposed  to  yield  a  solution 
f,  out  of  all  possible  solutions  in  solution  space  F,  which  gives 

a  minimvim  norm  |  |  (ij)  -  g  f)  |  |  .  It  seems  contradictory  to  say  that 
with  regularization,  another  solution  in  space  F  can  be  found 
which  gives  an  even  lower  value  of  the  norm.  It  is  noted  also 
that  in  all  cases  using  the  newly  defined  norm  (Equation  18)  the 
f ’  solution  contains  only  slight  deviations  from  the  f . 

Further  investigation  of  the  computer  routine  showed  the 
integration  form  assumed  in  Equations  (11a  and  b)  is  different 
from  the  form  used  in  the  back  calculation  (Equation  17) .  The 
N  values  at  the  low  a  values  not  only  contain  the  intrinsic  sum 
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of  least  squares  of  the  problem,  but  also  a  small  but  discernible 
part  due  to  the  difference  in  the  integration  forms  assumed  in 
the  forward  and  backward  mapping.  The  back  calculation  integra¬ 
tion  uses  the  Simpson's  Rule,  while  the  forward  mapping  uses  the 
Trapezoidal  Rule.  As  the  value  of  a  increases,  the  part  in  N 
attributable  to  the  different  integration  forms  will  decrease. 
With  further  increase  in  a,  the  regularization  begins  to  over¬ 
smooth  the  solution  form,  such  that  N  increases  again,  producing 
an  artificial  minimum  in  N. 

In  scheme  II,  the  backward  mapping  integration  routine  was 
changed  to  the  same  form  as  that  of  the  forward  mapping. 

Figure  2  shows  the  N  values  as  a  function  of  a  in  schemes  I 
and  II.  In  scheme  I,  a  minimum  is  observed,  but  the  effect  of 
regularization  on  the  solution  form  at  this  value  of  a  is  min¬ 
imal  comparing  with  the  solution  form  obtained  at  a=0.  In 
scheme  II,  as  expected,  N  is  an  increasing  function  with 
increasing  a. 

Using  the  relaxation  spectrum  given  by  Tobolsky  and 
Catsiff  [Reference  14] ,  the  three  sets  of  data  (in  log  form)  were 
generated  using  Equations  (1-3) .  All  data  contain  up  to  eight 
significant  figures.  Using  scheme  II,  all  sets  gave  very  good 
approximations  of  the  original  relaxation  spectrum  (Figure  3) 
without  regularization.  The  number  of  significant  figures  in 
the  data  sets  were  gradually  reduced,  and  the  unregularized 
solution  forms  in  all  cases  gradually  show  discernible  noise 
levels.  Finally,  a  "white  spectrum"  of  noise,  is  added  to  the 
two  significant  figure  data  of  log  modulus.  The  noise  is 
scaled  as  such  that  the  data  will  contain  2%  noise.  The  unreg¬ 
ularized  solution  form  is  shown  in  Figure  4a.  This  clearly 
demonstrates  that  quadratic  programmingalone  using  the  reduced 
norm  will  give  good  approximations  of  the  relaxation  spectrum, 
providing  that  the  initial  data  contain  sufficient  degree  of 
accuracy. 
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Figure 


.  The  N  values  as  a  function  of  a.  The  left  is 
calculated  with  Scheme  I,  and  the  right  with 
Scheme  II. 
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The  regularization  functional  was  used  in  an  attempt  to 
remove  the  noise  on  the  solution  form  obtained  from  the  data  with 
the  2%  noise.  The  scaling  factor  a  which  controls  the  relative 
contribution  of  the  regularization  was  gradually  increased.  At 
log  a  =  -19  (Figure  4)  the  big  valleys  in  the  middle  of  the  t 
ranges  were  diminished,  but  the  high  values  at  the  low  x  range 
also  suffered  from  over-regularization.  Several  n-values 
(Equation  12)  were  tested  and  all  gave  similar  results.  This 
is  not  really  surprising  considering  the  way  the  regularization 
term  is  written.  The  equation  is  regularizing  the  linear  form 
of  the  solution,  so  the  "real"  slopes  in  the  high  magnitude 
region  are  larger  or  comparable  to  that  of  the  noise  in  the  low 
region.  Regularization  of  the  low  region  necessarily  has  to 
over-regular ize  the  high  magnitude  region.  One  obvious  alternative 
is  to  regularize  only  the  region  that  contains  noise,  but  then 
a-priori  knowledge  of  what  the  solution  form  should  be  is  required. 

Another  alternative  is  to  rewrite  the  regularization  term 
such  that  the  function  to  be  regularized  is  log  H(t) ,  such  that 
the  "real"  slope  values  are  reduced.  But  the  term  will  still 
have  to  be  expressed  in  the  linear  form,  so  that  it  can  be  added 
to  the  term  N  for  the  quadratic  programming  routine.  The  way 
to  accomplish  that  is  not  immediately  obvious.  But  even  if  this 
can  be  accomplished,  one  is  still  faced  with  the  problem  of  know¬ 
ing  what  a  value  will  best  remove  the  noises  in  the  solution 
form,  yet  not  over-regularize  it. 

The  apparent  success  of  the  regularization  routine  to  sub¬ 
stantially  improve  the  solution  form  of  the  molecular  weight 
distribution  (MWD)  determination  may  be  due  to  two  factors: 
tl)  the  problem  is  more  "ill-posed";  (2)  and  the  solution  f° 
takes  on  only  a  small  range  of  values.  Not  every  equation  in 
the  form  of  Equation  (4)  is  "ill-posed";  it  depends  on  the  nature 
of  the  kernel.  Obviously  if  the  kernel  is  a  delta  function,  the 
problem  would  be  extremely  well  posed.  For  the  molecular  weight 
distribution  determination  problem,  the  problem  is  so  ill-posed 
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that  the  unregularized  solution  is  dominated  by  noise;  so  much  so 
that  it  almost  resembles  a  random  pattern  [Reference  11,12].  It 
is  possible  that  the  two  integration  forms  employed  in  the 
forward  and  backward  mapping  procedures  produce  a  much  bigger 
effect  on  the  back-calculated  value  of  N.  After  regularization 
has  svibstantially  improved  the  solution  form,  the  effect  from 
the  two  routines  is  minimized.  So  the  function  outputted  as 
f '  is  a  much  better  approximation  of  f“  than  f .  Also  in  the  MWD 
determination,  the  solution  is  in  units  of , mole-fractions ,  and 
the  numbers  do  not  vary  over  several  orders  of  magnitude  like 
the  relaxation  spectrum.  If  one  chooses  to  display  the  relaxation 
spectrum  in  a  linear  scale  rather  than  the  log  scale,  regularization 
will  appear  to  work,  too.  In  light  of  this  discussion,  it  seems 
interesting  to  reinvestigate  the  MWD  determination  with  the  two 
integration  routines  being  the  same,  but  unfortunately  that  is 
beyond  the  scope  of  the  present  work. 

Keeping  the  number  of  initial  data  points  (85)  constant,  the 
nvimber  of  spectrum  computed  points  was  varied.  With  the  2%  noise 
data,  the  noise  on  the  solution  form  decreases  with  decreasing 
number  of  spectrxim  points  (Figure  5)  ;  and  the  solution  obtained 
with  33  spectrum  points  became  a  reasonable  approximation  of  the 
original  relaxation  spectrum.  This  means  that  by  keeping  the 
output/input  data-point  ratio  low,  the  quadratic  programming 
can  compensate  for  the  noise  level  in  the  initial  input  data. 

This  is  verified  by  increasing  the  number  of  output  points  with 
the  set  of  data  without  noise.  The  acceptable  solution  form 
obtained  previously  became  noisy. 

Figure  6b  shows  the  solution  form  obtained  by  scheme  II 
using  Tobolsky's  experimental  data.  The  number  of  spectrum  points 
were  decreased  to  31  points  (with  85  points  of  input  modulus  data) . 
Figure  6a  shows  the  solution  form  obtained  with  scheme  I  under 
identical  conditions.  Obviously,  by  decreasing  the  number  of 
output  points  both  schemes  will  decrease  the  noise  level  in  the 
solution  form  (compare  with  Figure  1) .  Because  scheme  II  is 
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weighting  each  data  point  equally  in  its  function  to  be  minimized, 
under  identical  conditions,  scheme  II  gives  better  solution  than 
scheme  I.  Or,  it  can  be  said  that  scheme  II  can  tolerate  higher 
degree  of  noise  in  the  input  data. 

Since  the  ratio  of  the  output/input  points  is  a  determining 
factor  in  obtaining  reasonable  solution  form,  it  would  be  inter¬ 
esting  for  future  works  to  see  if  better  solutions  can  be  obtained 
from  a  given  set  of  data  by  artificially  increasing  the  number  of 
input  points  through  interpolation. 

Because  of  the  similar  shape  of  the  kernels  between  the 
storage  and  stress-relaxation  modulus  equations,  the  two  sets  of 
data  behave  very  similarly  in  the  quadratic  programming  routine. 
For  the  loss  modulus  data,  the  kernel  does  not  overlap  such  a 
wide  region  of  t ,  that  for  a  given  noise  level ,  the  output/input 
ratio  needed  to  give  a  reasonable  relaxation  spectrum  is  higher 
than  a  comparable  storage  or  stress  relaxation  data.  Since 
experimentally,  storage  and  loss  moduli  are  usually  obtained  with 
the  same  measurement,  it  would  be  interesting  for  future  works 
to  combine  both  sets  of  data  in  the  quadratic  programming  routine 
to  see  if  that  would  have  the  same  effect  as  increasing  the  nxmiber 
of  input  data  points. 
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SECTION  IV 
CONCLUSION 

The  regularization  quadratic  programming  approach  to  infer 
relaxation  spectra  from  experimental  mechanical  data  was  studied. 
It  was  foiind  that  the  expression  -  K(f^))  will  not  give  a 

minimvim  as  the  regularization  weighting  parameter  a  is  increased. 
A  modification  was  made  to  the  function  to  be  minimized  so  that 
all  data  points  will  be  weighted  equally  in  the  program.  Quad¬ 
ratic  programming  is  found  to  be  sufficient  to  infer  relaxation 
spectra  if  the  input  data  have  a  high  degree  of  accuracy.  When 
the  data  are  not  sufficiently  accurate,  regularization  will 
over-regularize  the  high  value  region  before  it  can  improve  on 
the  low  value  region  because  of  the  nature  of  the  relaxation 
spectrum.  Lowering  the  number  of  output  points  can  compensate 
for  the  noise  level  in  the  input  data,  and  will  allow  inferring 
relaxation  spectra  from  experimental  data. 
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